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Experimental  evidence  has  shown  that  composites  comprised  Si  and  Sn  nanoparticles  embedded  inside 
a  matrix  are  the  most  promising  next  generation  anodes  for  Li-ion  batteries.  This  is  due  to  the  ability  of 
the  matrix  material  to  constrain/buffer  the  up  to  300%  volume  expansion  that  Sn  and  Si  undergo  upon 
the  formation  of  lithium  rich  alloys.  Damage  still  occurs  at  the  nanoparticle/matrix  interface,  and  hence 
further  materials  design  is  required  in  order  to  commercialize  such  anodes.  Initial  theoretical  works  have 
predicted  that  low  volume  fractions  and  high  aspect  ratios  of  the  nanoparticles  result  in  a  greater  mechan¬ 
ical  stability  and  hence  better  capacity  retention.  The  most  important  design  parameters,  however,  such 
as  particle  size  and  spacing  have  not  been  considered  theoretically.  In  the  present  study,  therefore,  a 
gradient  enhanced  damage  model  will  be  employed  to  predict  that  damage  during  Li-insertion,  is  negli¬ 
gible  when  the  particle  size  is  20  nm,  and  the  interparticle  half-spacing  greater  then  1 .5  times  the  particle 
diameter.  Furthermore,  from  the  matrix  materials  considered  herein  graphene  is  predicted  to  be  the  most 
promising  matrix,  which  is  consistent  with  recent  experimental  data. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

In  trying  to  develop  Sn  and  Si  based  anodes,  for  secondary  Li- 
ion  batteries,  it  has  been  shown  that  configuration  and  geometric 
parameters  drastically  affect  their  electrochemical  performance. 
This  is  due  to  the  fact  that,  upon  the  formation  of  Li-rich  alloys, 
Sn  and  Si  expand  300%,  which  leads  to  fracture;  hence  the  effi¬ 
ciency  of  the  anode  depends  on  its  ability  to  withstand  these 
volume  changes.  Experimental  research  has  shown  that  fracture  is 
minimized  when  the  active  material  with  respect  to  Li  has  dimen¬ 
sions  in  the  nanoscale  [1,2],  while  further  mechanical  stability  is 
obtained  by  embedding  the  nanoparticles  in  a  matrix  in  order  to 
constrain  and  buffer  their  expansion  [3].  It  has  been  shown  for 
example  that  embedding  Si  nanoparticles  in  carbon  allowed  for 
a  capacity  retention  of  lOOOmAhg-1  for  twenty  cycles  [4],  while 
attaching  Si  nanoparticles  on  cellulose  gave  a  capacity  retention 
of  1400 mAhg-1  for  fifty  cycles  [5].  Furthermore,  embedding  Sn 
nanoparticles  in  carbon  allowed  for  a  capacity  of  500  mAhg-1  for 
100  cycles  and  attaching  Sn  on  cellulose  resulted  in  a  capacity  of 
500mAhg_1  for  forty  cycles  [6].  Not  all  Si  and  Sn  based  nanocom¬ 
posites,  however,  exhibit  such  good  capacity  retentions,  since  the 
volume  expansions  of  the  Si  within  the  composite  can  result  in 
severe  fracture  in  most  cases  [7,8].  As  a  result  the  capacity  decreases 
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significantly,  since  fracture  limits  the  connectivity  in  the  anode,  and 
active  material  can  be  lost  in  the  electrolyte.  Although  models  have 
been  developed  that  correlate  the  active  site  volume  change  to  the 
resulting  capacity  [9],  it  is  not  possible  to  obtain  a  material  selection 
that  provides  a  promising  capacity  without  a  significant  volume 
expansion  [10].  It  is,  therefore,  important  to  develop  theoretical 
models  that  can  predict  damage  in  nanocomposite  anodes. 

Linear  elastic  fracture  mechanics  have  been  employed  to  pre¬ 
dict  the  stable  configurations  [11]  and  to  develop  design  criteria  for 
selecting  the  matrix  material  and  determining  the  preferable  vol¬ 
ume  fractions  of  the  active  site  that  will  limit  fracture  [12].  These 
studies  indicated  that  increasing  the  volume  fraction  and  decreas¬ 
ing  the  aspect  ratio  of  the  Sn  or  Si  will  result  in  greater  instabilities 
during  fracture,  and  this  was  consistent  with  experimental  obser¬ 
vations  [13].  A  more  detailed  correspondence,  however,  between 
damage  and  capacity  showed  that  the  capacity  strongly  depends 
on  the  area  fraction  and  volume  average  of  the  active  particles 
[14].  Hence,  the  most  important  geometric  parameters  that  must  be 
considered  are  the  explicit  size  of  the  Sn  particles  and  their  spacing. 

In  doing  so,  standard  continuum  mechanics  studies,  as  those 
employed  thus  far  [11-13]  for  studying  anodes,  cannot  be  used, 
as  they  cannot  account  for  the  specimen  microstructure  and  scale. 
Instead  gradient  theories  must  be  utilized.  The  unique  feature  of 
gradient  theories  is  that  they  explicitly  consider  the  scale  of  the 
microstructure,  by  introducing  a  characteristic  length  unique  to 
the  material  at  hand.  Since  the  damage  that  occurs  in  Sn/C  anodes 
resembles  crumbling  [14],  i.e.  a  mode  of  distributed  damage,  the 
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most  appropriate  model  to  capture  this  behavior  is  that  of  gradi¬ 
ent  enhanced  damage  mechanics.  In  this  model  the  properties  of 
the  matrix  govern  the  damage  behavior  and,  therefore,  a  variety  of 
materials  will  be  considered,  similarly  as  in  [12],  in  order  to  deter¬ 
mine  which  material  selection  exhibits  minimum  damage.  Then  the 
extent  of  damage  will  be  computed  for  different  radii  of  the  Si  or 
Sn  particles,  providing  information  for  the  optimum  inter-particle 
spacing.  It  is  noted  that  in  fabricating  anodes  the  above  mentioned 
Sn  and  Si  based  nanocomposites  are  mixed  with  C  black  and  PVDF; 
in  the  present  study  the  effects  of  these  additives  to  the  mechanical 
stability  are  not  considered,  as  the  main  factor  resulting  in  capacity 
decay  is  fracture  of  the  Sn  or  Si  [11,14]. 


2.  Theoretical  approach 

2.1.  Damage  material  model 

One  defines  damage  as  the  alteration  of  any  material  physical 
property  due  to  the  presence  or  the  nucleation/growth  of  defects 
(microcracks,  voids,  delamination,  etc.).  In  order  to  characterize, 
represent  and  model  at  the  macroscopic  scale  the  effects  of  such 
distributed  defects  and  their  growth  on  the  material  behavior,  vari¬ 
ous  continuum  damage  models  have  been  developed.  These  models 
use  a  set  of  continuous  damage  variables  that  are  supposed  to 
describe  aspects  of  the  internal  material  structure  associated  with 
the  irreversible  (dissipative)  effects.  According  to  the  pretty  wide 
and  abstract  definition  of  the  phenomena  they  describe,  there  is  a 
variety  of  damage  variables  used  in  literature,  ranging  from  scalars 
over  first-,  second-,  forth-,  up  to  eighth-order  tensors.  A  survey 
on  this  subject  can  be  found  in  [15-18].  The  evolution  of  dam¬ 
age,  meaning  (i)  the  nucleation  of  new  microcracks  resulting  in 
distributed  microcracking,  as  well  as  (ii)  the  propagation  (growth) 
of  already  existing  microcracks,  induces  anisotropy  even  in  initially 
isotropic  materials.  However,  in  a  lot  of  problems  it  is  sufficient,  at 
least  from  the  phenomenological  point  of  view,  to  assume  that  the 
development  of  damage  does  not  affect  material  isotropy,  since  the 
microcracks  can  be  randomly  oriented  so  that  no  preferable  direc¬ 
tion  exists  or  the  microvoids  have  a  spherical  shape  and  isotropic 
distribution.  In  these  cases  scalar  damage  variables  adequately 
describe  the  local  state  of  a  damaged  material.  Two  classical  inter¬ 
pretations  of  this  approach  can  be  given.  The  first  one  (cf.  [15]) 
represents  the  ratio  between  the  area  dAD  of  the  intersection  of 
all  microcracks  and  microvoids  with  the  total  area  dA  of  the  plane 
section  D  =  dAD/dA.  The  second  one  (cf.  [16])  represents  the  current 
volume  fraction  of  the  voids  in  the  representative  volume  element 
D  =  dVPldV. 

From  both  definitions  it  follows  that  D  e  [0, 1  ],  where  D  =  0  stands 
for  the  undamaged  material  and  D  =  1  represents  maximum  dam¬ 
age  and  hence  complete  loss  of  integrity.  For  further  details  on  the 
subject  of  continuum  damage  mechanics  one  is  referred  to  [17,18]. 


gradient  parameter  which  implicitly  introduces  a  damage  internal 
material  length,  which  can  be  shown  to  be 


(2) 


and  H(p(d)  represents  an  appropriate  interaction  function.  The 
parameter  Ld  describes  the  size  of  the  damage  localization  zone, 
cf.  [19,23,24].  The  scalar  variable  d  measures  the  degree  of  mate¬ 
rial  stiffness  loss  and  is  monotonically  related  to  the  damage  ratio 
(reduction  of  the  Young’s  modulus) 


(3) 


where  Ec  stands  for  the  current  effective  Young’s  modulus  (or  better 
named  secant  stiffness  modulus)  that  corresponds  to  the  particular 
damage  state  and  E0  represents  its  initial  value.  Strictly  speaking, 
the  initial  material  always  contains  some  defects,  but  it  is  assumed 
that  these  are  accounted  for  in  the  virgin  material  properties.  The 
total  energy  n  is  obtained  in  standard  manner  as  the  difference 
between  the  potentials  of  internal  and  external  body  forces  b  and 
traction  forces  t  as 

n  =  njnt  -  next  =  /  d ,  <?d ,  v<pd)dv  -  f  u  (pb)dv 

Jq  Jo, 


u  t dA. 


(4) 


Field  equations  for  u  and  (pd  can  now  be  obtained  by  minimizing  n 
with  respect  to  these  variables: 


find  {u,  (pd }  =  argmin  {fl(u,  (pd) |u  =  u*  on  9£2U},  (5) 


where  the  internal  variable  d  is  considered  constant  with  respect 
to  variation,  for  details  see  [20-22]. 

Following  common  thermodynamic  considerations,  the  damage 
driving  force  is  defined  by  the  enhanced  free  energy  function  Eq. 
(1 )  as  the  conjugate  quantity  of  the  governing  variable: 


(6) 


The  evolution  of  the  damage  in  time  is  determined  by  a  damage 
potential  (pd  as 


d  =  ic  pL; 

Odd 


k  >  0,  <  0,  k(j)d  =  0, 


(7) 


where  d  represents  the  rate  of  the  damage  variable  and  k  stands 
for  the  rate  of  the  consistency  parameter.  The  detailed  form  of  (pd 
employed  in  this  work  is  given  in  Appendix  A.  Differential-algebraic 
inequalities  as  given  by  Eq.  (7)  are  common  in  the  modeling  of 
inealstic  materials,  see  e.g.  [25]. 

Eqs.  (5)  and  (7)  now  completely  determine  our  problem.  For 
details  on  the  solution  procedure,  see  once  again  [20-22]. 


2.2.  A  gradient-enhanced  damage  material  model 

In  order  to  account  for  the  scale  of  the  microstructure  [19],  a 
damage  material  model  including  a  regularization  strategy  based 
on  gradient  enhancement  [19]  of  the  free  energy  function  is  used 
in  the  present  study.  The  enhancement  is  formulated  through  an 
interaction  potential  [20-22],  which  depends  on  a  newly  intro¬ 
duced  variable  field  and  its  gradients 

(€:C:e)+g(d)+cJ-\\V<pci\\2  +  t!f[<pd-H<p(d)]2,  (1) 

where  pd  represents  the  energy  penalizing  the  difference  between 
the  corresponding  non-local  and  local  fields,  cd  stands  for  the 


2.3.  Geometries  considered 

Based  on  experimental  evidence  [4],  the  active  sites  are  treated 
as  spherical  inclusions  embedded  periodically  inside  a  matrix 
(this  approach  was  also  followed  in  previous  theoretical  modeling 
[28,11,3]),  as  illustrated  in  Fig.  1.  Considering  spherical  symme¬ 
try  allows  for  the  3D  problem  to  be  reduced  to  a  two-dimensional 
axisymmetric  configuration,  cf.  [29].  In  order  to  reduce  the  com¬ 
putational  effort  necessary  to  perform  the  damage  analysis,  the 
originally  rectangular  unit-cell  is  replaced  by  a  cylindrical  one,  see 
Fig.  1.  That  way,  a  small  error  is  introduced  which,  however,  is  cer¬ 
tainly  smaller  than  that  one  due  to  the  difference  between  regular 
and  random  distributions  and  thus  has  no  significant  influence  on 
the  results  in  general  [30].  The  unit  cell  problem  is  then  solved  for 
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periodic  boundary  conditions  ,  by  which  one  displacement  direc¬ 
tion  is  equal  to  zero  at  the  outer  boundary  of  the  unit  cell  as  shown 
in  Fig.  2a;  therefore,  the  effect  of  a  neighboring  nanoparticle  can 
be  accounted  for.  For  purposes  of  comparison  the  case  of  clamped 
boundary  conditions  (Fig.  2b)  is  also  considered  by  letting  both  dis¬ 
placement  components  be  zero  at  the  outer  boundary  of  the  unit 
cell.  The  geometry  of  the  representative  volume  element  and  con¬ 
sequently  the  matrix  geometry  is  assumed  to  be  cylindrical.  The 
volume  expansion  of  the  Li  inclusion  is  applied  as  a  uniform  dis¬ 
placement  in  the  radial  direction;  hence  a  300%  volume  expansion 
corresponds  to  a  58.74%  increase  of  the  inclusion  radius.  Therefore, 
it  will  be  assumed  that  the  diameter  increases  by  1/3  during  maxi¬ 
mum  Li-insertion,  which  corresponds  to  a  237.5%  volume  increase 
in  the  radial  direction.  The  behavior  of  a  material  at  such  large 
volume  changes  cannot  be  very  accurately  described  using  a  small- 
strain  material  model.  However,  it  can  give  first  insights  into  the 
problem  of  the  evolution  of  damage  in  the  matrix  material  sur¬ 
rounding  the  Li-active  inclusion.  The  analysis  formulated  in  the 
present  paper  is  based  on  large-strain  considerations. 

As  mentioned  in  Section  2.2  the  damage  model  employed  here  is 
sensitive  to  the  material  parameters  of  the  matrix,  hence  by  com¬ 
puting  the  damage  that  results  for  various  matrix  materials,  the 
most  promising  matrix  can  be  selected.  In  order  to  be  able  to  com¬ 
pare  with  the  predictions  of  earlier  studies  [12],  the  same  materials 
as  those  considered  in  [12]  (A1203,  B4C,  BeO,  WC,  [31  ])  will  also  be 
considered  here;  it  should  be  noted  that  such  materials  are  not 
used  in  actual  anodes.  Therefore,  in  addition  to  those  materials,  Cu 
[32]  and  graphene  [33]  will  be  treated  here,  since  they  have  been 
explicitly  employed  as  matrix  materials  in  nanocomposite  anodes. 
Cu  has  been  used  as  the  matrix  with  Sb  [35]  or  Sn  as  the  active 
sites  [34]  (upon  lithiation,  the  Sn-Cu  or  Sb-Cu  intermetallic  com¬ 
pounds  decompose  into  Sn  and  Sb  nanoparticles  dispersed  in  the 
inert  Cu  matrix).  More  recently,  graphene  has  been  recently 
employed  as  the  matrix  in  Sn-based  anodes  and  was  shown  to 


a _ tk-O _  b  u=0 


Fig.  2.  Numerical  analysis  model  including  periodic  (a)  and  clamped  (b)  boundary 
conditions. 


Table  1 

Material  parameters  for  various  ceramic  matrices  (E:  Young’s  modulus,  v:  Poisson’s 
ratio,  ft\  uniaxial  tensile  strength). 


Material 

E  (GPa) 

V 

ft  (MPa) 

A1203 

345.0 

0.23 

255.0 

B4C 

450.0 

0.21 

155.0 

BeO 

400.0 

0.24 

246.0 

WC 

700.0 

0.24 

345.0 

Cu 

130.0 

0.355 

210.0 

Graphene 

1000.0 

0.4 

130000.0 

Table  2 

Material  parameters  of  the  gradient  damage  model. 

cd  (GPa  nm2 )  /3d  (GPa) 

OL\ 

400.0  100.0 

1.5 

be  very  effective,  allowing  for  good  capacity  retentions  above 
550mAhg-1  [36-38];  it  is  noted  that  graphene  on  its  own  has  a 
low  capacity  300mAhg-1.  Although  Sn/graphene  anodes  are  not 
nanocomposites  in  the  traditional  sense,  as  graphene  has  a  lay¬ 
ered  structure,  based  on  TEM  images  [36]  it  can  be  seen  that  Sn 
nanoparticles  are  fully  surrounded  by  graphene  and  therefore  the 
configuration  employed  in  the  present  study  (nanoparticles  sur¬ 
rounded  by  a  matrix)  can  be  used  to  model  damage  in  such  systems. 
Furthermore,  the  authors  of  [36]  also  refer  to  such  Sn/graphene  sys¬ 
tems  as  nanocomposites.  An  experimental  measure  for  the  elastic 
modulus  of  graphene  does  not  exist  and  therefore  the  value  com¬ 
puted  by  MD  simulations  [33]  will  be  employed.  Inserting  these 
values  into  Eq.  (6)  provides  an  internal  length  of  2  nm. 

3.  Results  and  discussion 

To  illustrate  the  development  of  damage  for  the  various  mate¬ 
rial  systems  considered,  the  results  will  be  presented  by  plotting 
the  damage  ratio,  given  in  Eq.  (3),  over  the  radial  distance  from  the 
inclusion.  The  results  presented  are  obtained  for  the  volume  expan¬ 
sion  of  the  active  site  of  237.5%,  which  corresponds  to  an  increase 
in  diameter  of  the  inclusion  of  1/3  of  its  original  size. 

In  order  to  investigate  the  behavior  of  different  materials  which 
can  be  used  as  an  anode  matrix,  a  series  of  simulations  were  per¬ 
formed  by  fixing  the  radius  of  the  spherical  inclusion  at  1 00  nm  and 
the  half-distance  between  inclusion  centers  at  lOOOnm.  The  elas¬ 
tic  material  parameters  and  the  corresponding  tensile  strength  are 
listed  in  Table  1.  The  parameters  specific  to  the  damage  model  and 


Damage  distribution  for  different  materials-clamped  boundary-betad=100 


Fig.  3.  Damage  distribution  along  the  radial  distance  from  the  inclusion  center  for 
different  materials  of  the  surrounding  matrix. 


346 


B.J.  Dimitrijevic  et  al.  /  Journal  of  Power  Sources  206  (2012)  343-348 


used  in  the  gradient-enhancement  strategy  are  listed  in  Table  2; 
inserting  these  values  in  Eq.  (2),  gives  an  internal  length  of  2  nm. 

In  Fig.  3  it  is  seen  that  the  damage  is  always  highest  at  the  active 
site-matrix  interface  since  that  is  where  the  highest  tensile  stresses 
are  present,  and  as  the  distance  from  this  interface  increases  the 
damage  decreases.  From  Fig.  3  it  is  predicted  that  from  the  various 
matrix  materials  examined,  the  use  of  Cu  results  in  the  highest  dam¬ 
age,  while  graphene  in  the  least  damage,  which  would  be  expected 
considering  graphene’s  layered  structure,  and  the  resulting  ability 
to  accommodate  the  expansion  of  the  active  sites  (such  as  Si,  Sn). 
These  theoretical  predictions  are  in  agreement  with  experimental 
data,  since  when  Cu  is  used  as  the  matrix  in  Sn-Cu  anodes  the  capac¬ 
ity  retention  is  not  improved  over  the  pure  Sn  anodes,  indicating 
that  Cu  is  not  effective  in  buffering  the  Sn  volume  expansions  [34]. 
Using  graphene,  however,  as  the  matrix  in  Sn-C  anodes  allows  for 
stable  capacity  retention  over  500mAhg-1  [38],  suggesting  a  high 
mechanical  stability. 

Concerning  the  remaining  materials  considered  (A1203,  B4C, 
BeO,  WC),  it  is  interesting  to  note  that  B4C  is  the  most  promising, 
and  this  is  consistent  with  the  predictions  in  [  1 2  ].  It  should  be  noted 
that  in  the  model  employed  here  damage  depends  on  the  ratio  of 
the  material  parameters 

/t2(  1-V)  fR1 

1  £(l-2v)(l  +  v)‘  w 

Hence,  any  material  which  has  a  ratio  similar  to  that  of  graphene 
would  be  a  promising  matrix  candidate,  whereas  materials  with  r\ 
similar  to  that  of  Cu  should  be  avoided.  Of  course,  detailed  exper¬ 
imental  studies  should  be  performed  in  order  to  determine  the 
material  parameters  for  materials  at  the  nanoscale. 

Since  graphene  is  the  most  promising  matrix,  its  parameters 
will  be  used  in  the  sequel  as  we  proceed  to  examine  the  influence 
of  microstructure.  In  order  to  investigate  the  influence  that  particle 
spacing  has  on  anode  stability  a  new  set  of  simulations  is  performed 
in  which  the  radius  of  the  inclusion  is  kept  constant  at  100  nm,  but 
the  matrix  radius/half  spacing  between  neighboring  inclusions  is 
allowed  to  attain  the  values  of  200  nm,  300  nm,  500  nm  or  1 000  nm. 

In  Fig.  4  we  illustrate  that  the  damage  tends  to  0  once  the 
distance  from  the  inclusion  center  reaches  250  nm,  provided  the 
interparticle  spacing  is  large  enough.  However,  for  a  small  matrix 
radius  of  200  nm  damage  extends  to  the  next  neighboring  particle, 
and  therefore  the  anode  fractures  completely.  This  means  that  for 
nanoparticles  of  the  size  considered  here  a  larger  interparticle  dis¬ 
tance  gives  greater  mechanical  stability,  which  is  consistent  with 
the  earlier  prediction  that  lower  active  site  volume  fractions  allow 
for  higher  stability  [12].  Furthermore,  in  Fig.  4  it  is  seen  that  peri¬ 
odic  boundary  conditions  result  in  higher  damage  than  the  clamped 
boundary  conditions.  However,  it  can  be  concluded  that,  at  least  for 
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Distance  from  the  inclusion  center  [nm] 

Fig.  4.  Damage  distribution  along  the  radial  distance  from  the  inclusion  center  for 
different  half-distance  between  inclusion  centers  and  inclusion  radius  of  100  nm. 


this  model,  boundary  conditions  play  a  role  only  if  the  interparti¬ 
cle  half-distance  (matrix  radius)  is  quite  small  since  at  spacings 
larger  than  1.5  up  to  3  times  the  particle  diameter  both  boundary 
conditions  give  overlapping  results. 

In  order  to  better  illustrate  the  effect  of  interparticle  spacing,  the 
distribution  of  the  “material  damage”,  defined  as  percentage  of  the 
material  stiffness  loss  MD  [%]  =  (1.0  -fid))  x  100  for  two  different 
cases  of  interparticle  spacing  are  given  in  Fig.  5.  It  is  seen  that  the 
size  and  shape  of  the  damage  zone  is  highly  influenced  by  the  inter¬ 
particle  spacing  and  the  interaction  between  the  corresponding 
damage  zones.  Particularly,  for  small  interparticle  spacing,  when 
the  damage  zones  between  the  particles  interact  (according  to 
Fig.  4),  damage  is  highly  non-uniform,  whereas  for  large  interpar¬ 
ticle  spacing  damage  appears  to  be  radially  symmetric. 

In  addition  to  considering  particle  size  and  spacing,  it  is  of  inter¬ 
est  to  examine  volume  fraction  effects.  In  [12]  it  was  shown  that 
smaller  volume  fractions  result  in  a  greater  difficulty  for  cracks 
to  propagate.  In  Fig.  6,  however,  it  is  illustrated  that  keeping  the 
volume  fraction  constant,  while  varying  the  inclusion  radii  and  dis¬ 
tances,  results  in  significant  differences.  It  is  seen  that  reducing 
the  inclusion  radius  below  20  nm  causes  a  decrease  of  the  maxi¬ 
mum  damage,  but  proportionally  wider  spreading  of  the  damaged 
zone.  For  example,  although  the  damage  ratio  is  smallest  when  the 
inclusion  radius  is  20  nm,  50%  of  the  matrix  is  damaged,  whereas 
when  the  inclusion  radius  is  lOOnm,  the  maximum  damage  ratio 
is  higher,  but  only  25%  of  the  matrix  is  damaged.  This  phenomenon 
is  due  to  the  fact  that  in  the  former  case  the  interparticle  spacing 
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Fig.  5.  Damage  distribution  for  inclusion  radius  of  lOOnm  and  half-distance  between  inclusion  centers  of  133  nm  (a)  and  400 nm  (b)  for  the  case  with  periodic  boundary 
conditions. 
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Table  3 

Summary  of  results  when  graphene  is  used  as  matrix. 


347 


Particle  diameter  (nm) 

Inter-particle  spacing  (nm) 

Width  of  damage  zone  in 
matrix  starting  from  the 
interface  with  active  site  (nm) 

Damage  ratio  (initial  and  end 
values) 

20 

40 

10 

0.48-0.478 

20 

60-100 

20-25 

0.065-0 

40 

80 

20 

0.49-0.48 

80 

160 

40 

0.5-0.45 

200 

400-2000 

240 

0.6-0 

Distance  from  the  inclusion-matrix  boundary  [nm] 

Fig.  6.  Damage  distribution  along  the  radial  distance  from  the  inclusion-matrix 
boundary  for  different  inclusion  radii  and  distances,  and  fixed  volume 
matrix-inclusion  ratio  (first  number:  inclusion  radius,  second  number:  half¬ 
distance  between  inclusion  centers). 

becomes  comparable  to  the  internal  material  length  given  in  Eq.  (2) 
which  defines  the  size  of  the  damaged  zone. 

In  concluding  these  calculations  it  is  of  interest  to  produce  a  fig¬ 
ure  similar  to  Fig.  4,  but  for  a  smaller  particle  size.  In  Fig.  7,  therefore, 
the  inclusion  radius  is  kept  constant  at  10  nm  and  the  half-distance 
between  inclusion  centers  is  varied  between  20  and  50  nm,  to  allow 
for  similar  volume  fractions  as  in  Fig.  4.  In  comparing  Figs.  4  and 
7  it  is  seen  that  when  the  inclusion  radius  is  for  example  lOOnm 
the  maximum  damage  ratio  has  the  same  value  regardless  of  inter¬ 
particle  spacing,  whereas  when  the  inclusion  radius  is  lOnm  the 
maximum  damage  ratio  becomes  practically  negligible  once  the 
particle  half-spacing  is  1.5  times  larger  than  the  particle  diame¬ 
ter.  This  explicitly  illustrates  the  significant  effect  that  inclusion 
size  has  in  the  mechanical  and  hence  electrochemical  stability  of 
anodes.  Table  3  summarizes  the  predictions  that  the  theoretical 


model  gives  (in  Figs.  5-7)  for  a  graphene  matrix  with  different  par¬ 
ticle  sizes  and  interparticle  spacings.  It  is  noted  that  in  the  figures 
presented  the  particle  radius  is  stated,  along  with  the  radius  of  the 
surrounding  matrix,  while  in  this  table  the  spacing  between  two 
particles  is  given,  which  is  twice  the  matrix  radius. 

4.  Conclusions 

The  present  study  employs  a  gradient  dependent  damage  model 
to  explicitly  account  for  the  effect  that  Sn  or  Si  volume  fractions, 
particle  size  and  spacing  have  in  the  stability  of  nanocomposite 
anodes.  It  is  illustrated  that  when  the  interparticle  spacing  is  com¬ 
parable  to  the  particle  size  the  damage  is  non-uniform,  whereas 
when  the  interparticle  spacing  is  large  the  damage  is  localized  and 
radially  symmetric.  Furthermore,  it  is  shown  that  when  the  particle 
size  is  20  nm  and  the  interparticle  half-spacing  greater  than  30  nm, 
the  damage  during  Li-insertion  is  negligible.  Therefore,  it  can  be 
concluded  that  an  optimal  microstructure  would  consist  of  parti¬ 
cles  with  diameters  of  20  nm  or  below  and  interparticle  spacings 
at  least  1.5  times  their  diameter.  Finally,  it  should  be  noted  that  by 
considering  various  matrix  materials,  it  is  shown  that  least  damage 
occurred  when  the  ratio  of  the  material  parameters  specified  in  Eq. 
(8)  is  equal  to  that  obtained  for  graphene,  while  it  is  worst  when  it 
is  equal  to  that  for  Cu. 
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Appendix  A. 
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In  order  to  solve  the  differential-algebraic  system  Eq.  (7)  the 
evolution  equation  for  damage  is  discretized  in  time  employing 
a  Backward-Euler  scheme.  In  the  present  contribution  a  damage 
material  model  sensitive  to  tension  is  used.  To  achieve  this  task,  a 
decomposition  of  the  strain  tensor  into  a  “positive”  (e+)  and  “neg¬ 
ative”  (e-)  part  is  introduced 

3  3 

€+  =  ^2  ^  +  |e'^  ®  €  =  “  |e'^  ^ 
1  =  1  1  =  1 

where  £2-  stands  for  the  eigenvalues  of  the  strain  tensor,  while  fi2 
stand  for  the  corresponding  eigenvectors.  A  thermodynamically 
consistent  formulation  of  the  threshold  condition  dependent  on 
the  positive  part  of  the  strain  tensor  has  been  developed  by  [26] 
and  subsequently  used  by  e.g.  [27,21,23].  Here  it  assumes  the  form 


Fig.  7.  Damage  distribution  along  the  radial  distance  from  the  inclusion-matrix 
boundary  for  different  values  of  the  half-distance  between  inclusion  centers  and 
inclusion  radius  lOnm. 


4>i  :=  rn  -  2 


:  C  :  €~  +  2e+  :  C  :  €“)  <  0. 


(10) 
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In  the  numerical  examples  in  the  rest  of  the  paper  the  following 
softening  function  and  damage  potentials  are  used 

/(d)  =  (l-d)2,  g(d)=2-r,— H„(d) 

=  g(0)-g(d),  (11) 

where  r\  is  a  model  material  parameter  that  represents  the  damage 
threshold,  and  the  parameter  a\  controls  the  rate  of  softening. 
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